I hate this:(

Could we not have w = cholesky(X).U and work with a w::UpperTriangular instead of a dense matrix? Tried it locally, no real gain in performance though.

Could we not have w = cholesky(X).U and work with a w::UpperTriangular instead of a dense matrix?

Yeah, but this won't work with some of the AD backends (I know, it's super-annoying...). If you have a look at our compat-code for ReverseDiff (I believe), I think you'll see that we have to do some custom stuff to compute the pullback.

Tried it locally, no real gain in performance though.

I don't think we'd expect it to because internally we're iterating over the relevant elements of the matrix anyways, i.e. we're not gaining anything by telling the rest of the computational path that we're actually working on a lower-triangular matrix because it already assumes the given matrix is lower-triangular.

Completely unrelated, but if it is not 100%-guaranteed that you always end up with an upper triangular matrix when calling cholesky (which I think you can't if AbstractMatrix is supported), it would be better to work with .UL instead of .U (as we already do in other places of Bijectors and other libraries).

It's 100% guaranteed that it's available though, right? So it's a question of efficiency, not correctness.

The problem here is that

1. We need to work with a lower-triangular matrix (unless we re-implement _link_chol_lkj to work with vector).
2. We drop the type-information in the construction of w, hence we don't actually know if it was a uppper- or lower-triangular (in the case where we do cholesky(...).UL).

All in all, it seems we need an additional diversion, e.g.

link_chol_lkj(x::LowerTriangular) = link_chol_lkj(lower_triangular(parent(x)))
link_chol_lkj(x::UpperTriangular) = link_chol_lkj(transpose(upper_triangular(parent(x))))
link_chol_lkj(x::AbstractMatrix) = _link_chol_lkj(x)  # assume it's lower-triangular

?

I would guess it's mainly for efficiency reasons. But it's difficult to say if there are other implications as well, e.g., regarding AD. An alternative to .UL would be something like what's used in PDMats: https://github.com/JuliaStats/PDMats.jl/blob/fff131e11e23403931a42f5bfb3384f0d2b114c9/src/chol.jl#L6-L11 That should also be quite efficient and you could continue working with upper-triangular matrices.

I would guess it's mainly for efficiency reasons.

https://github.com/JuliaLang/julia/blob/db7971f49912d1abba703345ca6eb43249607f32/stdlib/LinearAlgebra/src/cholesky.jl#L515-L527

But it's difficult to say if there are other implications as well, e.g., regarding AD.

Hmm fair, though IMO something like this seems like it would be a bug with the AD package, no?

An alternative to .UL would be something like what's used in PDMats: https://github.com/JuliaStats/PDMats.jl/blob/fff131e11e23403931a42f5bfb3384f0d2b114c9/src/chol.jl#L6-L11 That should also be quite efficient and you could continue working with upper-triangular matrices.

Me and @harisorgn were just having a chat and we're thinking of replacing upper_triangular(parent(cholesky(X).U)) with

cholesky_upper(x) = upper_triangular(parent(cholesky(X).U))
cholesky_lower(x) = lower_triangular(parent(cholesky(X).L))

to make it less likely that we forget or mess up somewhere.

But we can make it

cholesky_upper(x) = upper_triangular(parent(PDMats.chol_upper(cholesky(X))))
cholesky_lower(x) = lower_triangular(parent(PDMats.chol_lower(cholesky(X))))

But are you sure there's not a good reason for why the default is copy? Of course it's more mem-intensive, but will stuff like LowerTriangular(U') lead to slower computation paths (since you're now working with adjoint(U) rather than something that is actually lower-triangular)? E.g. indexing adjoin(U) surely involves more computations than indexing copy(adjoint(U)).

The InverseFunctions.test_inverse that is called in here fails on some calls, but passes on others, even with the same sample x. CI on previous commits of this PR failed because of it. Seems weird, not sure why.

@harisorgn do we generate a random number somewhere inside the tests? If so, maybe fix the random seed?

There was no RNG between these lines, it was after sampling from LKJ. I believe the issue was numerical as there was no explicit zero-filling in the old Matrix-version link function, which caused mismatches. I've added zero-filling and tested it a bunch locally, should be fine 🤞